module exercise_lecture_1
use my_own_da
    implicit none
    PRIVATE derivative,Vexpt

	type vectorfield
	 type(my_taylor) F(2)
	end type Vectorfield

  INTERFACE OPERATOR (.d.)
     MODULE PROCEDURE derivative    
  END INTERFACE

  INTERFACE exp
     MODULE PROCEDURE Vexpt    
  END INTERFACE


contains 




  FUNCTION asint( S1 )
    implicit none
    TYPE (my_taylor) asint,tr,ti,ur,ui,ar,ai,temp
    TYPE (my_taylor), INTENT (IN) :: S1
	real(dp) tr0,ti0,xr,xi
    integer i

Tr=Sqrt(1.0_dp-s1**2)
tr0=Tr%A(0)
Tr%A(0)=0.0_dp
   
Ti=s1
ti0=Ti%A(0)
Ti%A(0)=0.0_dp

xr=tr0/(tr0**2+ti0**2)
xi=-ti0/(tr0**2+ti0**2)
   
ur=xr*tr-xi*ti
ui=xi*tr+xr*ti

ar=ur
ai=ui

tr=ur
ti=ui


   do i=2,3
    temp=ur*tr-ui*ti
    ti=ui*tr+ur*ti 
	tr=temp
     tr=-tr;ti=-ti;
	 
	 ar=ar+tr/i
	 ai=ai+ti/i
    
   enddo

   asint=asin(s1%a(0))+ai

  END FUNCTION asint


  FUNCTION asint_int( S1 )
    implicit none
    TYPE (my_taylor) asint_int,x,t,a
    TYPE (my_taylor), INTENT (IN) :: S1
	real(dp) tr0,ti0,xr,xi
    integer i,no,j(2)

     x=s1%a(0)+(1.0_dp.mono.1)

T=1.0_dp/Sqrt(1.0_dp-x**2)
     
	no=3
	a=0.0_dp
	x=(1.0_dp.mono.1)
	do i=0,no-1
     j=0;j(1)=i
	 a=a+ (t.sub.j)*x**(i+1)/(i+1)
	enddo

t=s1
t%a(0)=0.0_dp

    asint_int=asin(s1%a(0))
  	do i=1,no
     j=0;j(1)=i
	 asint_int=asint_int+ (a.sub.j)*t**i
    enddo


  END FUNCTION asint_int

  FUNCTION asint_inv_vecfield( S1 )
    implicit none
    TYPE (my_taylor) asint_inv_vecfield,t,ti,temp
    TYPE (my_taylor), INTENT (IN) :: S1
    TYPE (vectorfield) vec
	real(dp) cos0,sin0,asin0
    integer i,no,k,j(2)
    
asin0=asin(s1%A(0))
cos0=cos(asin0)
sin0=s1%A(0)

ti=(sin0*cos((1.0_dp.mono.1)/cos0)+cos0*sin((1.0_dp.mono.1)/cos0)-sin0)
temp=ti
t=(1.0_dp.mono.1)/cos0

    
	no=3

     vec%f(2)=0.d0;     
    do i=1,no
           vec%f(1)=-(ti-(1.0_dp.mono.1)) +vec%f(1) 
           ti=exp(vec,temp)
	enddo
    
           ti=exp(vec,t)

    asint_inv_vecfield=asin0
	temp=s1-s1%a(0)
	
	do k=1,no
	 j=0;j(1)=k
	 asint_inv_vecfield=asint_inv_vecfield+(ti.sub.j)*temp**k
    enddo



  END FUNCTION asint_inv_vecfield 


  FUNCTION derivative( T, I )
    implicit none
    TYPE (my_taylor) derivative
    TYPE (my_taylor) , INTENT(IN) :: T
    INTEGER, INTENT(IN) :: I
    INTEGER K1,K2,J(2),I1,I2

     derivative=0.0_DP
      I1=0
	  I2=0

	  IF(I==1) THEN
	    I1=1
	  ELSEIF(I==2) THEN
	  I2=1
      ELSE
       RETURN
	  ENDIF 

     DO K1=I1,3
	 DO K2=I2,3
     IF(K1+K2>3) CYCLE
      J(1)=K1;J(2)=K2;

      IF(I==1) THEN
	  derivative=(T.SUB.J)*DBLE(J(1))*(1.D0.MONO.1)**(J(1)-1)*(1.D0.MONO.2)**(J(2))+derivative
      ELSE
	  derivative=(T.SUB.J)*DBLE(J(2))*(1.D0.MONO.1)**(J(1))*(1.D0.MONO.2)**(J(2)-1)+derivative
	  ENDIF

	 ENDDO
	 ENDDO
     

  END FUNCTION derivative

  FUNCTION Vexpt( V,S1 )
    implicit none
    TYPE (my_taylor) Vexpt,T
    TYPE (vectorfield), INTENT (IN) :: V
    TYPE (my_taylor), INTENT (IN) :: S1
    integer i
    REAL(DP) NORMA,NORMB,EPS
    EPS=1.E-5

   T=S1
   Vexpt=S1
   NORMB=1.D30
   
   do i=1,100000
    
    T=(V%F(1)*(T.D.1)+V%F(2)*(T.D.2))/I
    Vexpt=Vexpt+T
    NORMA=FULL_ABS(T)
	IF(NORMA<EPS) THEN
	 IF(NORMA>=NORMB) THEN
	  RETURN
	 ELSE
	  NORMB=NORMA
	 ENDIF
	ENDIF
   enddo
    
  END FUNCTION Vexpt 

  FUNCTION FULL_ABS( S1 )
    implicit none
    REAL (DP) FULL_ABS
    TYPE (my_taylor), INTENT (IN) :: S1
    integer i

   
   FULL_ABS=0.D0

   do i=0,9
        FULL_ABS=FULL_ABS+ABS(S1%A(I))
   enddo
    
  END FUNCTION FULL_ABS 






  FUNCTION asint_inv( S1 )
    implicit none
    TYPE (my_taylor) asint_inv,t,ti,temp
    TYPE (my_taylor), INTENT (IN) :: S1
    TYPE (vectorfield) vec
	real(dp) cos0,sin0,asin0
    integer i,no,l,k,j(2),it
    
asin0=asin(s1%A(0))
cos0=cos(asin0)
sin0=s1%A(0)

ti=(sin0*cos((1.0_dp.mono.1)/cos0)+cos0*sin((1.0_dp.mono.1)/cos0)-sin0)
temp=ti
t=(1.0_dp.mono.1)/cos0

    
	no=3
   
     ti=(1.0_dp.mono.1)-temp+(1.0_dp.mono.2)
	   
	do it=1,no
     temp=0.d0
    do k=0,no
    do l=0,no
	 j=0;j(1)=k;j(2)=l;
     if(k+l>no) cycle
     temp=temp+ (ti.sub.j)*ti**k*(1.0_dp.mono.2)**l
	enddo
	enddo
	 ti=temp
	enddo
    ti=ti/cos0


    

	temp=(s1-s1%a(0))
    call print(temp,6)

	
	do k=1,no
	 j=0;j(1)=k
	 asint_inv=asint_inv+(ti.sub.j)*temp**k
    enddo

    asint_inv=asin0+asint_inv


  END FUNCTION asint_inv 


end module exercise_lecture_1
